Computer simulations of liquid silica: equation of state and liquid-liquid phase 

transition 
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We conduct extensive molecular dynamics computer simulations of two models for liquid silica [the 
model of Woodcock, Angell and Cheeseman, J. Phys. Chem. 65, 1565 (1976); and that of van Beest, 
Kramer and van Santen, Phys. Rev. Lett. 64, 1955, (1990)] to determine their thermodynamic 
properties at low temperature T across a wide density range. We find for both models a wide range 
of states in which isochores of the potential energy U are a linear function of T^/^ ^ as recently 
proposed for simple liquids [Rosenfeld and P. Tarazona, Mol. Phys. 95, 141 (1998)]. We exploit 
this behavior to fit an accurate equation of state to our thermodynamic data. Extrapolation of 
this equation of state to low T predicts the occurrence of a liquid-liquid phase transition for both 
models. We conduct simulations in the region of the predicted phase transition, and confirm its 
existence by direct observation of phase separating droplets of atoms with distinct local density and 
coordination environments. 

PACS numbers: 65.50.+m 64.30.+t, 64.60.My, 64.70. Ja 



I. INTRODUCTION 

First order liquid-liquid phase separation, in which two 
liquids of distinct chemical composition coexist, are com- 
mon in multicomponent systems. However, there has in 
recent years been a growing interest in first order liquid- 
liquid phase transitions that occur without a change of 
composition, but rather with a change in density p as 
temperature T or pressure P is varied. Experimental 
evidence for the occurrence of such transitions has been 
found for a wide range of systems, including Si pi, 2|, 
I, Se, S J, |, AI2O3-Y2O3 melts |, C @, H2O § |, 
and P pf7 Liquid-liquid transitions have also been ob- 
served in molecular dynamics (MD) computer simula- 
tions of Si [|0[ |l|, H2O [H m and C (l|. The- 
oretical studies have long predicted liquid-liquid tran- 
sitions for a variety of model fluids; for examples, see 
Refs. @ [6|, ^, ^, ^ ^, H, g, g . 

In the case of water, the proposed liquid-liquid phase 
transition occurs in the supercooled liquid, i.e. for T 
less than that of the melting line . Closely associated 
with the possibility of a liquid-liquid transition in su- 
percooled water is the phenomenon of polyamorphism in 
the amorphous solid occurring below the glass transition 
temperature, Tg. Polyamorphism refers to the occurrence 
of distinct amorphous solid forms of a material p6| . 
In the mostprominent cases of polyamorphism, such as 
water [|7[ |2^, an abrupt first-order-like transition 
occurs from a low-density form to a distinct high-density 
form as the amorphous material is compressed at low T. 
For water it was proposed that the observed polyamor- 
phism of the amorphous solid is a sub-Tg manifestation of 
the thermodynamic instability associated with the liquid- 
liquid phase transition ]T2t . 

Computer simulation studies of the ST2 water 
model [3(J] support this view of the relationship between 



liquid-liquid phase transitions and polyamorphism. The 
qualitative features of water polyamorphism are clearly 
displayed in simulations of ST2 water |3l|. Correspond- 
ingly, the critical temperature Tc marking the onset of 
liquid- liquid phase separation has been determined [ p^ , 
as well as a characteristic pattern of thermodynamic 
"precursors" for T > Tc consistent with the instability at 
T = Tc . These precursors include the occurrence of a 
density maximum and a compressibility maximum. Sim- 
ulation studies of the thermodynamic properties of two 
other water-like models, TIP4P H and SPC/E Q, find 
the same pattern of thermodynamic precursors observed 
for ST2 for T > Tc H, H, and also display polyamor- 
phism in simulations of the amorphous solid. However, 
in these systems, simulations of the liquid at lower T, 
to test for the onset of a liquid-liquid phase transition, 
have not yet been attempted due to prohibitively long 
equilibration times. 

Substances that are structurally similar to water, such 
as Si, Ge, Ge02, and Si02 (silica), have the potential 
to exhibit similar behavior |3q] . Particular attention has 
been paid to silica because of its technological and geo- 
logical importance. Polyamorphism is indeed observed in 
compression experiments on amorphous silica |3^, ^ , 
and is also qualitatively reproduced in computer simula- 
tions [^0|, ^ . Though not as dramatic as is found for 
amorphous solid water, the polyamorphism of silica may 
also be due to a trend toward liquid-liquid phase sepa- 
ration |4^, Q. Indeed, liquid state simulations of the 
silica model of Woodcock, Angell and Cheeseman (de- 
noted here "WAC silica" ) |Q have shown that the same 
pattern of thermodynamic precursors of the liquid-liquid 
phase transition found in water simulations also occurs 
in this system @. However, as in TIP4P and SPC/E 
water, the low T simulations required to test for an ex- 
plicit liquid-liquid phase transition in WAC silica have 
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not been attempted to date. 

At the same time, recent advances concerning the prop- 
erties of Uquids at and below the melting line are expand- 
ing our ability to study the states where liquid-liquid 
phase transitions may occur. Of particular importance 
is the recent prediction |^ that an isochorc of the po- 
tential energy U should be a linear function of T^/s Jq^. 
a simple, cold, dense liquid. Since this relation is pro- 
posed to be valid in the limit of low T, it is a useful 
relation for studying the properties of a deeply super- 
cooled liquid. Notably, several recent works showed that 
this prediction is obeyed at low T for binary Lennard- 
Jones liquids This observation provided a 

physical basis for extrapolating U (and thermodynamic 
properties derived from it) to T near Tg, and so made 
possible determination of the Kauzmann temperature for 
this system. 

In this paper, we examine in detail the behavior of two 
simulation models of silica: (i) WAG silica; and (ii) the 
widely used potential of van Beest, Kramer and van San- 
ten iQ, denoted here as "BKS silica." Our goal is to 
determine if either model displays a liquid-liquid phase 
transition. We find that our computer simulation data 
for the thermodynamic properties of these silica mod- 
els obey the prediction of Ref. over a wide range of 
T and V. We exploit this result to construct equations 
of state for BKS and WAG silica, and find that liquid- 
liquid phase transitions are predicted for both models at 
low T. We then conduct simulations near the predicted 
phase transition, and confirm the occurrence of liquid- 
liquid phase transitions for both BKS and WAG silica by 
direct observation. 



II. MOLECULAR DYNAMICS SIMULATIONS 

The calculations presented here for WAG silica are 
based on the data set generated for Ref. |^ . These simu- 
lations consisted of iV = 450 atoms (300 O, 150 Si atoms) 
and were conducted in the constant {N, V, E) ensemble. 
{E is the total internal energy, V is the volume.) The ef- 
fects of electrostatic interactions were incorporated using 
the Ewald summation technique [ pT| . 

For the BKS model, we conduct new simulations of a 
system of iV = 1332 atoms. As for WAG silica, Ewald 
summations are used to include electrostatic interactions. 
The Ewald parameter (a in the notation of Ref. ||5l| ) is 
fixed to 2.5 nm~^ for all state points simulated. For each 
state point, the system is equilibrated to near the desired 
T using periodic velocity rescaling. All averages are re- 
ported for constant {N, V, E) simulations that follow the 
equilibration stage. In all cases, averages are evaluated 
over a time that is at least ten times longer than the 
average time required for an Si atom to diffuse 0.2 nm. 

The BKS potential has the unphysical feature that the 
interaction energy of a Si and O atom pair diverges to 
—oo as their separation goes to 0. Though not a prob- 
lem at ambient T and P, this feature will occasionally 



manifest itself at high T and P. We have added a short 
range term to the BKS potential that prevents this from 
occurring, but which does not alter the form of the BKS 
potential at larger separations |52| . 

The (y, T) coordinates of the state points simulated 
for this work are shown in Fig. |^. Table | gives V for 
each of the isochores studied. 



III. TEMPERATURE DEPENDENCE OF 
POTENTIAL ENERGY ISOCHORES 

Ref. predicts that the isochoric T dependence of 
the potential energy U oi a simple, dense, cold liquid is 
given by, 

t/ = a + 6T3/^ (1) 

where a and b are constants for a given V. 

Here we test if Eq. |l| is obeyed by WAG and BKS silica. 
We plot isochores of U against T^/^ and fit a straight 
line to the data (Fig. ||). At all V studied, we find that 
Eq. 1^ fits the data well. Gonsistent with the prediction of 
Ref. the best fits occur for the smallest V, and the 
quality of the fits decreases somewhat as V increases. For 
BKS silica, we find that Eq. |l| fits to all isochores within 
numerical uncertainty. For WAG silica, most of the data 
available to us fits Eq. ^ within numerical uncertainty. 
However, for the largest V, systematic deviations from 
Eq. occur if all WAG data up to the highest T are 
included. Yet even for these large V isochores, the WAG 
data are consistent with an approach to the behavior of 
Eq. |l| at low T. By excluding several of the highest T 
data points, shown as open circles in Fig. ^ we recover 
a fit within numerical uncertainty even for the largest V 
isochores for WAG silica. 

The V dependence of the fit parameters a and b so 
obtained for both models is shown in Fig. The success 
of the fits to Eq. |l| over a wide range of T and V, and 
the smooth variation of a and b with V show that the 
predictions of Ref. appear to be valid for silica, a 
"complex" liquid with anisotropic molecular interactions, 
at least in the limits of low T and low V. 

In addition, it is important to note the V dependence 
of a. a provides an estimate of U in the limit T ^ 0, 
which classically is coincident with the limit as T — > 
of A, the Helmholtz free energy. For both WAG and 
BKS silica we find a range of V in which the curvature 
of the a versus V curve is negative. In this range, the 
condition for thermodynamic stability for a single phase, 
{d'^A/d^V)T > 0, is not satisfied Q, suggesting that 
both WAG and BKS silica would undergo a liquid-liquid 
phase separation at low T, if not preempted by crystal- 
lization or vitrification. 
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IV. MODEL EQUATION OF STATE 



That is, 



Having identified a region of validity for Eq. |l|, we can 
use this relation to construct a representation of the ther- 
modynamic properties of WAC and BKS silica in terms of 
a continuous function of T and V . The model equations 
of state so generated will allow us to clarify the liquid- 
liquid phase separation suggested by the V dependence 
of a in Fig. ^. 

We first fit polynomials to the V dependence of aiV) ~ 
Hn=Q'^nV'^ and h{V) ~ Z]t=o/^"^"' obtain a func- 
tional representation of U : 



U{V,T) = a{V)+b{V) T^/^ 



(2) 



The coefficients a„ and /3„ are given in Table ||. The in- 
ternal energy is E = U+Uk', for the classical ionic models 
of silica considered here, the kinetic energy is Uk = ^RT, 
where R is the gas constant. Hence, our model for E is. 
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E{V, T) = a{V) + h{V) T^/^ -|- -RT. 



(3) 



Next, we seek a functional representation of the en- 
tropy S{y,T), so that a model for A{V,T) can be ob- 
tained from A — E ~ TS. S at arbitrary V and T, rela- 
tive to the entropy S{Vq, Tq) of a reference state, can be 
evaluated by thermodynamic integration j53|. We carry 
this out in two steps, first along an isotherm, and then 
along an isochore. 

We compute the change A5t = S{V,Ta) - 5(Fo,To) 
along an isotherm T = Tq using A A = A£' - TqASt, 
where AA = A{V, To) - A{Vo,To) and AE = E{V, To) - 
E{Vo,To). A£; is evaluated from Eq. |. AA is a difference 
due to a volume change from Vq to ^ at T = Tq, and is 
found from an isothermal integration of dA = —PdV; 
that is, 



AA- 



I P{V',To) 

JVo 



dV. 



(4) 



To obtain a functional representation for AA therefore 
requires one for P{V) at T = Tq. We obtain the re- 
quired data from our MD simulations and fit to them a 
polynomial P ~ X]n=o7"/^"' where the density p= 1/V 
(Fig. |). The coefficients 7„ are given in Table WR along 
with the choices of Tq for WAC and BKS silica. Using 
this polynomial model for P{V,To), the integration in 
Eq. U yields an expression for AA which combined with 
that for AE, gives a model function for ASt [Q- In 
terms of E and P, the expression for ASt is, 



ASj 



1 



EiV,To) ~ E{Vo,To) 



P{V',To)dV' 



Vo 



(5) 

The change ASy = S{V,T) - S(y,To) is the entropy 
difference at fixed V due to a temperature change from 
To to T. This we find from an isochoric integration of. 



dT. 



(6) 



This evaluation is carried out using our representation of 
E{V,T) in Eq. |. 

Combining the contributions of both isothermal and 
isochoric changes, S at an arbitrary state point is given 
by. 



S{V, T) = SiVo, To) + ASt + ASy- 



(8) 



Using the model functions for E and S, we thus ob- 
tain a function modeling A(V,T). The equation of state 
PiV,T) is fomid from. 



PiV,T) 



dA 
dV 



(9) 



Note that the resulting expression for P{V, T) does not 
contain the unknown reference entropy S(Vo,To) since 
this constant disappears after the differentiation in Eq. ||. 

To summarize, we construct a model P(V, T) equation 
of state using as input, polynomial fits of (i) the V de- 
pendence of a and 6, and (ii) one reference isotherm of 
P. As a check of this equation of state, we compare in 
Fig. II isochores of P versus T, evaluated directly from 
simulation, and as calculated from the above modeling 
procedure. 



V. THERMODYNAMIC BEHAVIOR OF WAC 
AND BKS SILICA 

For the description of the thermodynamic properties of 
tetrahedrally coordinated liquids such as silica or water, 
it is useful to determine the location and shape of curves 
in the space of P, V and T at which specific thermody- 
namic conditions are met. In the present context, three 
such curves are important: 

(i) Along the "temperature of maximum density" 
(TMD) line the condition. 



dP 
df 



0, 



(10) 



is satisfied [Q. At such a point, an isobar of p as a 
function of T is a maximum, and at lower T, p decreases 
as T decreases. The presence of a TMD line is a hallmark 
of liquids in which local tetrahedral order is prominent, 
and is observed experimentally in silica, as well as in 
water. 

(ii) The metastability limit of the liquid, or spinodal 
line, is defined by |36[, 



dP 
dV 



= 0. 



(11) 
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(iii) Along the "K™'^^ line," the isothermal compress- 
ibility Kx is a maximum with respect to V at constant 
T. It is found by locating points satisfying, 



dKr 



= 0, 



where, 



}_(dV_ 



(12) 



(13) 



and then checking to confirm that the extremum so iden- 
tified is a maximum p6| , |55[ |. 

Spinodal lines are necessarily associated with a second 
order critical point that terminates a line of first order 
phase transitions |Q. When a such a critical point oc- 
curs aX T — a, K^^^ line will emanate from the critical 
point for T > Tc. However, the occurrence of a K^'^^ line 
does not imply the occurrence of a critical point at lower 
T; i.e. the occurrence of a K^^^ line is a necessary, but 
not a sufficient condition, for the occurrence of a critical 
point ||||. 

The locations of these lines are shown in Fig. |l| pro- 
jected onto the V-T plane; and in Fig. ^, projected onto 
the T-P plane. We find that the pattern of behavior re- 
vealed is qualitatively the same both for WAC and BKS 
silica, and that this is the same pattern found also from 
water simulations employing the ST2, TIP4P and 
SPC/E models ||. 

Most significant is the occurrence both for WAC and 
BKS silica, of a spinodal line in the low T liquid regime 
that is distinct from the liquid-gas spinodal boundary. 
This spinodal is the metastability boundary associated 
with a liquid-liquid phase transition. As T — > it co- 
incides with the points of inflection in the a versus V 
curves of Fig. |^. The critical point of this liquid-liquid 
phase transition occurs at the point of maximum T on 
the spinodal line. 



VI. LIQUID-LIQUID PHASE SEPARATION 

If the spinodal curves predicted by the equations of 

are correct, and if we can 



state developed in Section IV 



conduct equilibrium simulations in the unstable regions 
so identified, we should observe characteristic signs of 
phase separation. When a liquid-liquid phase separation 
occurs in a constant- simulation such as employed here, 
both phases coexist within the simulation box, each in 
its own distinct region, separated by an interface. Each 
phase will have a distinct bulk density as well as a distinct 
local structure. 

We therefore carry out new MD runs for both WAC 
and BKS silica at state points that approach the spin- 
odal curve associated with the liquid-liquid phase tran- 
sition. These state points are identified as squares in 
Fig. I. To facilitate comparison of BKS and WAC silica, 
we simulate a system of = 750 atoms for both. The 



lowest T state simulated for each model (T = 2000 K for 
BKS, T = 4000 K for WAC) is near the predicted critical 
point. Because of the low T, we are are unable to bring 
these lowest T states fully into equilibrium; however, the 
results obtained do serve to establish the trend in the 
behavior at the lowest T. 

To test for the occurrence of two phases with distinct 
local structure, we examine the local coordination en- 
vironment of the silicon atoms. We consider g(r), the 
Si-Si radial distribution function (RDF); A'Kr'^g(r)dr is 
the probability that a Si atom will be found at a dis- 
tance between r and r -I- dr of a reference Si atom. We 
decompose g(r) according to the contributions made by 
successive nearest neighbors (nn) of a given Si atom, la- 
belled in ascending order of distance from an Si atom. 
That is, we define sub-RDF's gi{r) according to 



(14) 



where Airr^ gi{r)dr is the probability that the iih. nearest 
neighbor of a randomly selected Si atom will be found at 
a distance between r and r -\- dr. 

Fig ^ shows (?5(r) for BKS and WAC silica for sev- 
eral different T along the isochores indicated by open 
squares m Fig. 0. For T above T^, 55 (r) is a unimodal 
function of r. As T decreases, the width of the 35 (r) dis- 
tribution increases. That is, rather than finding a more 
sharply defined 5th nn coordination environment as T de- 
creases, the distribution of locations of 5th nn's becomes 
broader. For T near Tc, g^{r) becomes bimodal. This 
behavior shows that two distinct populations of 5th nn 
coordination environments are emerging in the liquid as 
T decreases. Similar behavior is observed for ge{r), and 
more weakly in 57 (r) and gsif)- This is in contrast to 
the behavior of gi{r) through 54 (r) which simply be- 
come sharper and narrower distributions with decreasing 
T (not shown). 

If the emergence of these distinct coordination envi- 
ronments corresponds to the onset of liquid-liquid phase 
separation, we should find that Si atoms with similar co- 
ordination environments are spatially correlated. That 
is, there should be relatively compact droplets of the two 
distinct phases. 

To test for this, we require systems of larger spatial 
extent than the N = 750 atom simulations used to calcu- 
late 35 (r) . We therefore initiate simulations of a system 
of A^ = 6000 atoms for the lowest T states where phase 
separation should be most prominent (T = 2000 K for 
BKS, T = 4000 K for WAC). As above, excessive equili- 
bration times prevent us from bringing these states fully 
into equilibrium, and so we only use these simulations 
to establish the trend in behavior. We note that an in- 
completely equilibrated system should underestimate the 
amount of phase separation, since droplets (if present) 
will have had less time to form and grow. 

We examine snapshots of these N = 6000 atom BKS 
and WAC configurations. The minimum occurring be- 
tween the two peaks of 35 (r) for the lowest T in Fig. |^ 
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provides a convenient threshold r* for partitioning the Si 
atoms into two populations according to their 5th nn co- 
ordination. We show in Fig. || snapshots of the positions 
of all the Si atoms in the BKS and WAC systems. Light 
spheres are Si atoms having their 5th nn at a distance 
greater than r*; these atoms have a low-density coordina- 
tion environment, compared to the average. Dark spheres 
are Si atoms having their 5th nn at a distance less than 
r*; these atoms have a high-density coordination envi- 
ronment. Spatially correlated droplets of atoms with the 
same type of coordination environment are readily vis- 
ible, for both WAC and BKS silica. This is consistent 
with the occurrence of liquid-liquid phase separation for 
both of these sihca models. 



VII. DISCUSSION 

The two principal conclusions of the present work are: 

(i) For two silica models, BKS and WAC, there exists a 
wide range of T and V within which isochores of U for the 
liquid phase conform to Eq. This occurs in spite of the 
fact that the prediction of Ref. was made for simple 
liquids, and not for liquids with non-trivial local structure 
such as occurs in silica. This result suggests that the 
physical basis for Eq. |^ is quite robust and provides a 
valuable tool for probing the low T properties of a wide 
range of liquid systems. 

(ii) The model equation of state constructed by ex- 
ploiting Eq. |l| predicts the occurrence of a liquid-liquid 
phase transition in BKS and WAC silica. We confirm 
the presence of this transition by direct simulations near 
the predicted critical point. Thus BKS and WAC silica 
join the rank of simulation models for tetrahedrally coor- 
dinated liquids in which a liquid-liquid phase transition 
has been directly observed. Whether or not such a phase 
transition occurs in real liquid silica (see below), it's pres- 
ence in the behavior of BKS silica is important to note, 
since this model is currently in wide use for simulation 
studies of silica under a variety of conditions. 

Given the common behavior found for BKS and WAC 
silica, it is appropriate to inquire whether we should 
therefore expect to find the same pattern of thermody- 
namic behavior, including a liquid-liquid phase transi- 
tion, in real silica; or whether these two models share a 
common flaw that makes them unrealistic in this respect. 
To attempt to address this, we seek a basis for compar- 
ing the behavior of the BKS and WAC models with each 
other and with other tetrahedral liquids, both simulated 
and real. To proceed we choose as a scaling temperature 
T* the highest value of T reached along the TMD line. 

As shown in Table [II , the ratio T^/T* is approximately 
0.4 for both WAC and BKS silica. Assuming that this 
ratio is also valid for real silica gives Tc = 730 K, a tem- 
perature well below Tg for silica at 1450 K. Hence, if the 
WAC and BKS models are representative of the thermo- 
dynamic properties of silica, then we should not expect to 
directly observe an equilibrium liquid-liquid phase tran- 



sition in supercooled liquid silica. 

This is in contrast to the case of water. The ST2 model 
gives Tc/T* = 0.80, implying that = 194 K for real 
water. This is well above Tg = 136 K for water, though as 
yet still outside the easily accessible experimental range. 
(See however, Ref. ||^.) 

The present analysis is consistent with recent interpre- 
tations [^2[ ^ Q of the phenomenology of polyamor- 
phism as observed in amorphous solid silica and water. In 
this interpretation, the manifestation of polyamorphism 
for a system where the liquid phase exhibits the tendency 
toward a liquid-liquid phase transition at low T will de- 
pend on the relationship of Tg to Tc- When Tc > Tg, 
as may be the case for water, polyamorphic behavior 
in the amorphous solid will be prominent, with a large 
and sharp density increase observed when the glass in 
compressed. When Tc < Tg, behavior characteristic of 
polyamorphism will be weaker, though not necessarily 
absent, as is the case for silica. 

However, the precise nature of the behavior to be ex- 
pected when Tc < Tg remains an open area of research. 
For T < Tg, a standard application of equilibrium ther- 
modynamics is not appropriate and so predictions based 
on the model equations of state presented here do not ap- 
ply. Nonetheless, new thermodynamic approaches, based 
on a separation of the configurational degrees of free- 
dom (frozen in at Tg) and vibrational degrees of freedom 
(which are always in thermal equilibrium) may be em- 
ployed to derive a free energy expression that can be 
used used to locate the (cooling-rate dependent) location 
of the critical point In this approach, a negative 

curvature of the V dependence of U along the Kauzmann 
line would be sufhcient to guarantee the presence of an 
instability in the glassy phase. 

We also note that the value of T for the fC™"^ line at 
atmospheric P is about 80% of T* for both WAC and 
BKS silica (see Fig. ^). Applying this ratio to real silica 
suggests that a maximum of Kt (or, equivalently, a min- 
imum of the bulk modulus) could occur when liquid silica 
is compressed isothermally at T just above, but near Tg. 
This raises the possibility to directly observe a if™"^ line 
in experiments of supercooled liquid silica under pressure. 
Confirming or refuting this prediction would be an impor- 
tant step in establishing the applicability of the pattern 
of thermodynamic behavior presented here, to real silica 
and other tetrahedral liquids. 
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TABLE L Volume V and density p of each isochore simulated 
for BKS and WAC silica. The labels are used to identify 
isochores shown in the figures. 
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model 


n 


a„ 




7n 




WAC 





-9.13919 


-0.00627597 


122.772 


To 


= 7000 K 


1 


-38.8581 


0.108913 


-328.537 






2 


163.650 


-0.494217 


311.922 






3 


-306.599 


0.983621 


-328.537 






4 


205.422 


-0.683688 


29.1778 






5 


- 


- 


2.31542 




BKS 





-1.01756 


0.00144202 


479.401 


To 


= 5000 K 


1 


-61.1632 


0.0488104 


-708.648 






2 


293.988 


-0.381123 


409.488 






3 


-621.091 


1.069764 


-117.045 






4 


478.866 


-0.969089 


16.6264 






5 






-0.904994 



TABLE II: Coefficients of polynomial fits for V dependence 
of a and b; and the p dependence of P for T = To. The units 
for each coefficient are appropriate to give a in MJ/mol, b in 
MJ/(molT^/^), and P in GPa, when V is given in cm /g and 
p is measured in g/cm^. 





T* 


Tc 


T, 


TJT* 


Tc/T* 


Tc/T, 


WAC 
BKS 

ST2 


9000 
5000 
330 ill] 


4000 
2000 
235 @ 


<4000 
1380 [|| 
<235 


<0.44 
0.28 
<0.71 


0.44 
0.40 
0.71 


>1 
1.45 

>1 


SiOa 
H2O 


1823 59) 
277 |61] 


~ 730 
~ 194 


1450 [p0| 
136 |62] 


0.80 
0.49 


~ 0.40 
~ 0.70 


0.50 
1.4 



TABLE III: Comparison of characteristic temperatures of 
simulated and real tetrahedral liquids. T* for real Si02 and 
H2O are estimated as the ambient P values of the TMD; the 
actual value of T* for these systems is likely to be slightly 
higher. Upper bounds on Tg for WAC silica and ST2 water 
are taken simply as the lowest T at which equilibrium simula- 
tions of the liquid have been conducted. Values preceded by 
"~" are found by assuming that the values of Tc/T* found 
from simulations apply to the real substances. References are 
provided in brackets for those values not determined here. 




FIG. 1: Simulated state points and thermodynamic features 
of (a) BKS and (b) WAG silica. The model equations of state 
derived in Section IV are fit to MD simulation results obtained 



at the {V, T) points indicated by filled circles. In (b) open 
circles are WAG state points at which substantial deviations 
from Eq. |l| are observed, and so are excluded from the data 
set used to construct the equation of state. The BKS and 
WAG equations of state give estimates for the projection into 
the {V,T) plane of the spinodals (solid lines), TMD line (dot- 
dashed), and if™"^ line (dashed), as defined in Section ^ 
The open squares locate (V, T) points at which we test for 
liquid-liquid phase separation, as described in Section 




FIG. 2: Isochores of U versus T^^/^ for (a) BKS and (b) WAG 
silica. Symbols are U values obtained from MD simulation, 
while solid lines show linear fits to each isochore. To better 
view each isochore and their fitted lines, we plot U — kV, with 
k — 1 (g MJ)/(cm'' mol) and V measured in cm'^/g, so that 
each isochore is subject to a I'-dependent shift to separate it 
from the others. In (a) isochores Bl through BIO are shown 
from bottom to top; in (b) isochores Wl through W15 are 
shown from bottom to top. Also in (b) are shown high T 
points (open circles connected by dashed lines) on the largest 
V isochores which deviate from linear behavior and so are 
excluded from the fits. The statistical error for the U values 
shown in (a) and (b) does not exceed ±0.004 MJ/mol. 
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FIG. 3: V dependence of (a) a and (b) b for BKS and WAG 
silica, found from the linear fits to the isochores shown in 
Fig. 0. Solid curves are fits to the data of a fourth order 
polynomial in V . In (a) D is a scale factor to permit both 
curves to be compared in a single plot; for BKS, D = 1 and 
for WAG, D = 0.45. 
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FIG. 4: Isotherms of P versus V for BKS (To = 5000 K) and 
WAG (To = 7000 K) silica. Solid curves are fits to the data 
of a fifth order polynomial in p. 
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FIG. 5: Isochores of P versus T for (a) BKS and (b) WAG sil- 
ica. Symbols are values obtained from MD simulation, while 
lines are determined from the model equations of state devel- 
oped in Section For BKS, isochores Bl through B8 are 
shown from bottom to top; for WAG, W4 through W15 are 
shown, from bottom to top. 




FIG. 6: Estimates for the projection into the (T, P) plane of 
the spinodals (sohd lines), TMD hne (dot-dashed), and if""^ 
line (dashed), evaluated from our model equations of state for 
(a) BKS and (b) WAG silica. 
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FIG. 7: 55 (r) for (a) BKS and (b) WAG silica evaluated at 
the state points indicated by open squares in Fig. |l[ The 
system size used is A'' = 750 atoms. The arrows indicate 
the values of r* used to identify silicon atoms having high- 
density coordination or low-density coordination. For BKS 
silica, r* = 0.355 nm; for WAG silica, r* — 0.360 nm. 



(a) BKS 



(b) WAC 



FIG. 8: Snapshots of MD configurations of BKS silica (T = 
2000 K, p = 0.34 g/cm^) and WAC silica (T = 4000 K, p = 
0.41 g/cm^). The configurations are viewed face-on to one 
side of the simulation cell. There are N = 6000 atoms in 
each system, but only silicon atoms arc shown. Dark spheres 
are silicon atoms having high-density coordination; i.e. a 
5th nn silicon within a distance r < r*. Light spheres are 
silicon atoms having low-density coordination: i.e. a 5th nn 
silicon at a distance r > r*. The clustering of atoms with 
similar coordination into droplets is consistent with the onset 
of liquid-liquid phase separation. 



